% function [x,x2,x4] = model_x2(a,b,x0,tspan)
% Calculate <x>,<x^2> and <x^4> where x subject to a first order stochastic
% ODE:
%   \dot{x} = ax+ sigma*N
% 
% Analytical solution
% 
a = -1;
sigma = 1;
x0 = 1;

t = 0:0.1:5;
x = exp(a*t)*x0;
x2 = exp(2*a*t)*(x0^2+sigma^2/(2*a))-sigma^2/(2*a);
x4 = exp(4*a*t)*(x0^4+3*x0^2*sigma^2/a+3*sigma^4/(4*a^2))-exp(2*a*t)*(3*x0^2*sigma^2/a+3*sigma^4/(2*a^2))+3*sigma^4/(4*a^2);

close all;
figure(1);
plot(t,x,'LineWidth',2);
set(gca,'fontsize',20);
xlabel('time(s)');
ylabel('<x>');
print('-dpng','meanX');
print('-depsc2','meanX');

figure(2);
plot(t,x2,'LineWidth',2);
set(gca,'fontsize',20);
xlabel('time(s)');
ylabel('<x^2>');
print('-dpng','meanX2');
print('-depsc2','meanX2');

figure(3);
plot(t,x4,'LineWidth',2);
set(gca,'fontsize',20);
xlabel('time(s)');
ylabel('<x^4>');
print('-dpng','meanX4');
print('-depsc2','meanX4');

figure(4);
plot(t,x2,'r-',t,sqrt(x4-x2.^2),'b--','LineWidth',2);
set(gca,'fontsize',20);
xlabel('time(s)');
legend('<x^2>','std(x^2) = sqrt(<x^4>-<x^2>^2)');
print('-dpng','meanX2_varX2');
print('-depsc2','meanX2_varX2');